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Abstract 

We present theoretical and numerical studies on stiff, linear polyelec- 
trolytes within the framework of the cell model. We first review analytical 
results obtained on a mean-field Poisson-Boltzmann level, and then use 
molecular dynamics simulations to show, under which circumstances these 
fail quantitatively and qualitatively. For the hexagonally packed nematic 
phase of the polyelectrolytes we compute the osmotic coefficient as a func- 
tion of density. In the presence of multivalent counterions it can become 
negative, leading to effective attractions. We show that this results from a 
reduced contribution of the virial part to the pressure. We compute the os- 
motic coefficient and ionic distribution functions from Poisson-Boltzmann 
theory with and without a recently proposed correlation correction, and 
also simulation results for the case of poly(para-phenylene) and compare 
it to recently obtained experimental data on this stiff polyelectrolyte. We 
also investigate ion-ion correlations in the strong coupling regime, and 
compare them to predictions of the recently advocated Wigner crystal 
theories. 



^email: markusOchem.ucla. edu 
^email: holin@inpip-niainz.nipg.de 



1 Introduction 



"Polyclcctrolytcs arc polymers bearing ionizable groups, which, in polar sol- 
vents, can dissociate into charged polymer chains (macroions) and small coun- 
terions" [1] . The combination of macromolecular properties and long-range elec- 
trostatic interactions produces an impressive variety of phenomena [2, 3, 4]. It 
makes these systems interesting from a fundamental as well as a technological 
point of view. 

A thorough understanding of polyclcctrolytcs has become increasingly im- 
portant in biochemistry and molecular biology. This is due to the fact that 
virtually all proteins, as well as the DNA, are polyelectrolytes. Their interac- 
tions with each other and with the charged cell-membrane are still far from 
being understood. For instance, a puzzling question is why two equally charged 
objects should attract each other in the first place [5]. 

In this article we focus on stiff linear polyelectrolytes, which we subsequently 
approximate by charged cylinders. This is a relevant special case, applying to 
quite a few important polyelectrolytes, like DNA, actin filaments or micro- 
tubules. We treat the solvent in dielectric approximation and explicitely de- 
scribe only the small ions. Within Poisson-Boltzmann (PB) theory [6] and on 
the level of a cell model the cylindrical geometry can be treated exactly in the 
salt-free case [7, 8, 9, 10, 11, 12], providing for instance new insights into the 
phenomenon of the Manning condensation [13, 14]. 

The paper is organized as follows: First we review some essential features of 
the PB mean-field solution of the cell model. Then we discuss the applicability of 
PB theory to the ion distribution functions and show under which circumstances 
PB theory fails quantitatively (underestimated condensation) and qualitatively 
(overcharging, charge oscillations and attractive interactions). Next we present 
measurements of the osmotic pressure for the nematic phase of hexagonally 
packed polyelectrolytes as a function of density and compare it to PB predic- 
tions and the Manning limiting law. In the next section we study the particular 
case of poly(para-phenylene) by means of PB theory, including a correlation 
correction of the basis of a recently proposed Debye-Hiickel-Hole-cavity theory 
(DHHC) [15], and simulational results. The results are compared to recent ex- 
perimental data on this system [16, 17, 18]. We find that correlation effects 
enhance condensation and lower the osmotic pressure, yet are not fully able to 
explain the discrepancy to the experimental data. At the end we try to shed 
light onto the role of specific ionic correlations. Two-dimensional correlations 
on the surface of the rods are found to be present, but weakly developed. No 
hexatic order of the ions is observed. 

2 Simulation method and model system 
2.1 The Langevin thermostat 

We utilize molecular dynamics (MD) simulations using a Langevin thermostat 
[19] to study the equilibrium properties of our model system within the canonical 
ensemble. Technically this is achieved by integrating the stochastic differential 
equation 

= -VUi{r,})-~Tr,+Ut) (1) 
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on the computer, where rrii are the masses of particles with coordinates 
subject to a potential energy function U, a friction T linear in the velocity and 
a stochastic white noise ^j(t). The latter two can be thought of as imitating 
the presence of a surrounding viscous medium responsible for a drag force and 
random collisions, respectively. Since both effects have the same origin, they are 
related to each other by a simple version of the fluctuation-dissipation-thcorem. 
This can be exploited to choose their strength such as to converge towards the 
canonical state: 

(e,(t))=0 and {$M-ijit'))=QkBTTd,,S{t~t'), (2) 

where is the thermal energy. 

We use a standard Verlet integrator [20] to integrate (1). Time step St and 
friction coefficient T were set to 0.01 and 1.0 in Lennard-Jones units (see below). 

2.2 Interaction potentials 

We use the purely repulsive Lcnnard- Jones (LJ) potential to give the countcrions 
an excluded volume: 
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< r < rcut = 21/6^ 
: rcut < r. 



The advantage of including the —r~^ contribution instead of merely using the 
purely repulsive r~^^ is that Eqn. (3) is exactly zero beyond rcut and merges 
smoothly to this value at rcut, allowing a larger time step. 
The Coulomb potential of a charge Q is written as 

peoVcir) = (Q/eo)^, (4) 
r 

where £b = Pe^/Aire is the Bjcrrum length, the distance at which two unit 
charges have an interaction energy equal to 1//3 := fceT, eo is the positive 
elementary charge, and e is the product of the dielectric constants of vacuum, 
ffo, and the medium, Er, respectively. The Lennard- Jones unit a is always used 
as unit of length and e is always set to the thermal energy k^T. Temperature is 
implemented via the Bjerrum length. Mass is irrelevant and set equal to one - 
it would only be needed to translate the Lennard-Jones time tlj = a\Jm/e into 
"real" time. We intend to model an aqueous environment at room temperature, 
which implies £b ~ 7.14 A. 

Within the periodic boundary conditions employed during the simulations, 
the presence of such long-range interactions poses both mathematical and tech- 
nical difficulties. We use the most efficient FFT accelerated Ewald sum, the P^M 
code, which scales almost linearly with the number of charges [21, 22]. 



2.3 Generating a cell-geometry 

Compared to the spherical cell model, the cylindrical one presents one additional 
but crucial complication: the charged rod is infinitely long. Several methods have 
been proposed in the literature to handle this problem [23, 24, 25, 26]. They 
essentially all use as a unit cell a hexagonal prism with a certain height. This 
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approximates the cylindrical cell by a space-filling object. In the present work 
we take a cube of side length Lb and place a rod along the main diagonal. Upon 
periodically replicating this system the diagonal rod becomes infinitely long 
and an infinite triangular array of such replicated rods emerges. The resulting 
Wigner-Seitz cell of this lattice is a regular hexagon, which can alternatively be 
viewed as the unit cell. This is illustrated in Fig. 2. Observe that the symmetry 
of the replicated system is still cubic. The radius i? of a circle with the same 
area A is then given by i? = Lb/V7r\/3. This value is most appropriately used 
for comparing results between the hexagonal and the cylindrical cell model. 

If the line charge density of the rod is A, electroneutrality requires the number 
of w-valent countcrions to be = ^/SL^X/veo- Hence the average counterion 
density is given by n = N/L^ ~ ^/SX/vcqL'^ and is thus inversely proportional 
to instead of L^. The number of countcrions can therefore be written as 
N = {^/3X/veo)^^^n~^^'^ , which implies that a smaller density requires more 
particles. While this makes the simulation of dilute systems rather expensive, 
it gives at the same time rather small dense systems. The latter problem can 
be circumvented by combining blocks of 2x2x2, 3x3x3 or even more 
elementary cubes to a big cube and using the latter as the unit box for the 
periodic boundary conditions. 

The ratio between £d = (47r^B'y^7T,)~^/^ (the average Debye length of the 
countcrions) and the rod separation diod can be written as 



where the definition of the Manning parameter ^ Afe/eo was used, see also 
Sec. 3. Obviously, for strongly charged rods the Debye length is smaller than 
the separation of the two rods. Even for = 1 it is only half as large as the 
distance between rod axis and Wigner-Seitz boundary, and the neighboring cells 
effectively decouple. This statement is independent of density. 

3 Poisson-Boltzmann Essentials 
3.1 The Poisson-Boltzmann Solution 

At this stage we want to briefly recall the necessary knowledge about the PB 
solution of the cylindrical cell model [7, 9, 10, 11, 12]. Consider an infinitely 
long cylinder of radius ro and line charge density A > 0, coaxially enclosed in 
a cylindrical cell of radius R. Global charge neutrality of the system is ensured 
by adding an appropriate amount of oppositely charged w-valent countcrions. In 
the following only the case of no extra salt will be discussed. 

Within PB theory the individual countcrions are replaced by a cylindrical 
counterion density n(r), where r denotes the radial distance from the cylin- 
der axis. Defining the reduced electrostatic potential y and and the screening 
constant k > as 
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the PB equation can be written as 
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It will be useful to define the dimensionless Manning parameter ^ = A^b/so 
which counts the number of unit charges on the rod per Bjerrum length. In the 
following the main focus will be on the strongly charged case characterized by 
> 1. Fixing the two boundary conditions for the electric field, 

y'iro)^-2^v/ro and y\R) = 0, (8) 

the correctly normalized solution to Eqns. (7,8) can be written as [7, 11, 12] 

y{r) = -2 In I + 7"" cos (7 In | . (9) 

Insertion of the general solution from Eqn. (9) into the boundary conditions from 
Eqn. (8) yields two coupled transcendental equations for the two integration 
constants 7 and Rm- 

7 In — — = arctan — and 7 In = arctan — . (10) 

Rm 7 7 

Subtracting the left part from the right part in Eqn. (10) eliminates Rm and 
provides a single equation 

7 In — = arctan h arctan , (11) 

ro 7 7 

from which 7 can be obtained numerically. The second integration constant Rm, 
which will be referred to as the Manning radius, can be obtained from either 
one of the Eqns. in (10) as soon as 7 is known. Note finally that k and 7 are 
connected via 

K2i?2 ^ 2(1+7^). (12) 

Since the ^ and v only enter in the combination ^w, changing valence or 
electrostatic interaction strength is equivalent on Poisson-Boltzmann level. In 
particular, at given cell geometry {ro,R} the Manning radius Rm is a unique 
function of ^i;. Eqn. (11) implies the sequence of inequalities ln(i?/ro) < n/j < 
lTi{R/rQ)+^v/{(,v—l). The resulting asymptotic for 7 in the dilute limit i? — > 00 
gives rise to what is often called the "Manning limiting laws" . 

The integrated probability distribution of finding a mobile ion at distance 
r within a cylinder of radius re [rg ; R] can be determined analytically by 
integrating the charge density: 



P(r) 



1 f'^ 1 V 

- dr'27rr'eovn{r') = (l-—)+^tSin(^\n^). (13) 



This is the fraction of countcrions found within a cylinder of radius r. P{r) will 
be referred to as counterion distribution function. At r = Rm the last term in 
P vanishes, giving the Manning fraction /m '■— 1 — l/^w of ions within Rm- 
These are the ions which can not be diluted away, if the cylinder radius is sent 
to infinity. This phenomenon is referred to as Manning condensation [13, 14]. 
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3.2 How to quantify counterion condensation 

If the counterion distribution function P is known, the condensed counterion 
fraction can be characterized in the foUowing "geometric" way: Eq. (13) shows 
that P viewed as a function of ln(r) is merely a shifted tangent-function with 
its center of symmetry at (ln(i?M); /m)- Since, however, tan"(0) = 0, Manning 
radius and Manning fraction can be found by plotting P as a function of ln(r) 
and localizing the point of inflection. 

This property of P, derived within the framework of PB theory, can in turn 
be used to define the condensed fraction [27, 12]. It provides a suitable way to 
quantify counterion condensation beyond the scope of PB theory, and it is exact 
in the salt-free PB limit. From here on this method will be referred to as the 
inflection point criterion. It has the advantages oi (i) not fixing by definition the 
amount of condensed counterions (/m and i?M can be determined independently 
of each other), (ii) reproducing the salt-free PB limit, and (Hi) quantifying the 
breakdown of the coexistence of condensed and uncondensed counterions in the 
high salt limit. 

4 Generic ion distribution functions 

On the plain PB level the radius rg of the rod is not a completely independent 
variable, since it enters only in the combination R/rQ, see Eqn. (11). The rod 
in our simulation is built up from a sequence of small Lennard- Jones particles 
lined up along the main diagonal of the simulation cell. The distance of closest 
approach of two Lennard-Jones particles can then be identified with the radius 
To of the rod and yields an effective rod radius of a. 

Upon leaving the PB level, the ratio between counterion diameter a and 
rod radius rp becomes relevant, which has for instance been investigated in 
Ref. [25, 28]. Here we do not intend to perform a systematic investigation of 
such effects, but instead present later in Sec. 6 results for systems mapped to 
physically relevant parameters, DNA and two kinds of poly(p-phenylenes), in 
which the rod has a considerably larger diameter than the counterions. 

4.1 Density dependence within monovalent systems 

At fixed rod radius of rg the relevant variable R/rQ is changed by varying the 
cell size R and therefore the density. Here we present results for such a density 
scan for systems with €b/c € {2,3}. The monovalent and positively charged 
particles forming the rod were placed along the main diagonal with a center- 
center distance of 1.04245(7, giving a line charge density of A = 0.959279 60/(7. 
In connection with the two presented values for the Bjerrum length this yields 
Manning parameters of ^ = 1.919, 2.878 and corresponding Manning fractions of 
/m = 1 — 1/C = 0.4788, 0.6525 in the monovalent case. The cell radius R has been 
varied between 2.06(7 and 124(7, which for monovalent counterions corresponds 
to average ion densities of 7.2 x lO^^a^^ and 2.0 x lO^^cr^^, respectively. The 
inflection point criterion from Sec. 3.2 has been used to determine the radial 
extension of the condensed layer, and the fraction of ions within it, = 
P(i?^), where the star denotes the measured quantities. A graphical illustration 
of the distribution functions is given in Fig. 3. 
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The important things to observe are the fohowing. The measured condensed 
fraction is always larger than the fraction predicted by the PB theory. How- 
ever, the fraction from the simulation decreases monotonically with decreasing 
density towards the Manning limit /m- Such a deviation is to be expected, since 
PB is essentially a low density theory [29, 30]. 

In contrast to the clear tendency of the measured condensed fraction to de- 
crease upon dilution, the behavior of the condensation radius appears to be 
more complicated. There does not seem to exist a simple monotonic convergence 
of towards Rm- Rather, for high densities the measured condensation dis- 
tance is larger than the Manning radius, while for the investigated low densities 
it is smaller. Unfortunately, a clear-cut statement is difficult since the localiza- 
tion of the point of inflection in P as a function of In(r') is only possible with 
an error estimated to be of the order of 1%. 

For Bjerrum length £b = 1 cr, corresponding to ^ = 0.9593 < 1, the radial 
distribution functions are compared to the PB predictions in Fig. 4. It can be 
seen that the measured and PB predicted distributions almost coincide for all 
cell sizes under investigation. While this is to be expected for low densities, the 
remarkable agreement at high densities is somewhat surprising. Apparently, PB 
theory is a fairly good description of weakly charged systems, i.e., ones which 
are below the Manning threshold. 

4.2 Fitting to a generalized Poisson-Boltzmann distribu- 
tion 

Since the measured Manning radius R^ can be smaller than the PB prediction, 
but > /m and R^ < Rm is incompatible with PB theory, it is impossible to 
describe such a measured P{r) with a PB distribution function in which a some- 
what enlarged effective Manning parameter ^* is used. Nevertheless, it would 
be desirable to use at least some of the knowledge about the ion distribution 
function from PB theory for evaluating relevant observables in a simulation 
or in a real experiment. To the lowest order, e.g., within a two state model, 
the increase of counterion condensation can be ascribed to such an effective 
= 1/(1 — /^)- To break the monotonic connection between /m and Rm one 
assumes that the functional form of P is given by the PB form in Eqn. (13), 
but neglects the relation in Eqn. (11). Using Eqn. (10?), this suggests fitting the 
measured distribution in a region around the inflection point to the form 

It* / r l-£*v\ 
Pm{r) = 1 - — + tan UMn — + arctan ^ (14) 

with the three fit parameters 7* and rj.^ The condensed fraction is then given 
by = 1 — l/£,*v and the condensed radius by R^^ = exp{arctan[(^*w — 
l)/7*]/7*}. Fig. 5 illustrates this procedure for one system. While this approach 
might seem to be very powerful, it suffers from the drawback that the actual 
numbers depend on the chosen fitting region. Nevertheless, it provides an inde- 
pendent way of quantifying condensation and can even be applied to very noisy 

■^It might seem strange to regard the rod radius as a free parameter, but there are two 
reasons for this. First, this does not force the fit to coincide with the PB form outside the 
chosen fitting-region. Second, not even in a real experiment is the radius of a macroion always 
known in advance. Rather, it is often necessary to obtain this information within the same 
measurement [16]. 
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data, for which the locahzation of an inflection point in P as a function of hi(r) 
is otherwise virtually impossible. 



4.3 Multivalent ions 

The universal dependence of the ion distribution on the product is an artifact 
of Poisson Boltzmann theory. Fig. 6 shows examples of systems with different 
Manning parameters but the same value of ^v. Not only is the condensation 
enhanced as compared to PB theory, but the enhancement is stronger for the 
case involving multivalent ions. For instance, in the dense trivalent system the 
condensation is enhanced by 30% with respect to the PB prediction. 

The reason is that at given charge density multivalent systems have a lower 
number density and thus fewer particles. This lowers any kind of excluded vol- 
ume interactions, and also makes the entropic contribution less pronounced. 
Moreover, the relevant variable for describing the strength of correlation effects 
is £bv'^ [15]. Hence, an increase in valence is more important than an increase 
in Bjcrrum length. 

Within PB theory the expression for the contact density n{rQ) in the limit 
of infinite dilution is given by 



where eo^" is the surface charge density of the rod. For high rod charge this 
expression becomes independent of valence. Hence, replacing monovalent coun- 
terions by multivalent ones reduces their total number in the cell, but for a highly 
charged rod not their density at the rod surface. This result is a consequence 
of the contact value theorem and thus holds beyond the Poisson-Boltzmann 
approximation [9] . 

4.4 Addition of salt 

The influence of salt on the distribution functions has been discussed within the 
PB theory in detail in Ref. [12]. The general finding was that a low salt content 
leaves the picture of Manning condensation qualitatively unchanged, while at 
increasing salt concentration a crossover between Manning condensation and 
simple salt-screening occurs. 

If one starts adding N salt molecules to the salt free solution one finds 
that the inflection point gets shifted to smaller values of r, hence the layer of 
condensed counterions contracts. The second finding is that the amount of con- 
densation is only marginally increased. From a certain N on two more inflection 
points appear close to the cell boundary. This happens typically for a corre- 
sponding Debye length being of the order of the cell size itself, indicating the 
appearance of a characteristic, salt induced, change in the convexity of P (as 
a function of In r) . Upon a further increase in N one of the two new inflection 
points shifts towards smaller values of r, finally fusing with the Manning inflec- 
tion point and "annihilating" it. Roughly speaking, the inflection points vanish 
if the Debye length characterizing the salt content becomes smaller than the ra- 
dius of the condensed layer. In this case it is no longer meaningful to distinguish 
between condensed and uncondensed counterions. Indeed, for a very high salt 




(15) 
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content, where the Debye length is much smaller than the radius of the rod, the 
solution of the PB equation would be the one of a charged plane and one may 
consider all excess countcrions being condensed - no matter what the charge 
density of the rod is. 

For three systems with numbers of salt molecules TV £ {0, 104, 3070} a sim- 
ulation has been performed and compared to the PB prediction in Fig. 7. As in 
the salt-free case the computer simulations show a more pronounced condensa- 
tion effect towards the rod. Nevertheless, the shape of the distribution functions 
remains qualitatively the same. Note in particular that the appearance and dis- 
appearance of two points of inflection at = 104 and N = 3070 respectively, 
which leads to extremely small curvatures in the PB distribution functions, 
also leads to very straight regions in the measured distribution functions. The 
crossover from Manning condensation to screening, as described within the PB 
theory, can be expected to be essentially correct. It should not be overlooked, 
however, that the addition of salt ions will affect the distribution function much 
more dramatically if their valence is larger than the valence of the countcri- 
ons. In this case, the salt counterions will accumulate in the vicinity of the rod 
at the expense of the lower valent "real" countcrions as has previously been 
demonstrated in a system with a valence mixture of counterions [31]. 

The PB approach fails to describe the physical situation if one or more of the 
following conditions apply: (i) the electrostatic interactions are strong, (ii) the 
counterions are multivalent or (in) the density is high. Results of a simulation 
under such conditions can be seen in Figure 8. Here, a system with box length 
Lh/o' = 36.112, i.e., 60 monovalent counterions and a cell size of R/a ~ 15.481, 
and Ib/ct = 4.1698, i.e., ^ = 4, has been investigated after adding 1000 molecules 
of a 2:2 salt. Since this gives almost 17 times as much salt as counterions and 
a salt Debye length of £u/cr = 0.33 <C -R/cr, this can essentially be viewed as 
a charged rod in a bulk 2:2 electrolyte. The most characteristic feature of the 
charge distribution function is that it overshoots unity, showing a charge reversal 
of the rod at distances around r w 1.5 cr, while the simple PB prediction is clearly 
qualitatively off. This phenomenon is usually referred to as overcharging and has 
been predicted for the rod geometry first from hypernetted chain calculations 
[32] and later by a modified PB approach [33, 34] , and has also been observed in 
planar geometries [35]. Since P{R) = 1 for the reason of global electroneutrality, 
the overshooting above 1 at small distances implies the existence of a range 
of r- values at which the mobile ion system is locally positively charged, i.e., 
with the same charge as the rod, such that P{r) can eventually decay to 1. 
This is seen in the right frame of Fig. 8, which shows that n+2{r) > n_2(r) at 
r ^ 2a. Since P(1.5 a) « 1.45, the rod and its innermost layer of condensed ions 
could be viewed as an effective rod of radius 1.5 tr which is negatively charged 
with Manning parameter ^ = 1.8. Since this value is again larger than 1, it 
entails ion condensation, but this time of positive ions. In fact, it even leads to 
a second overcharging, as can clearly be seen in Fig. 8, where P{r) - in decaying 
from 1.45 - overshoots the value of 1 again. Overcharging can thus give rise 
to layering. In the presented example no less than three layers can clearly be 
made out. These local charge oscillations also reflect themselves in oscillations 
of the electrostatic potential, as demonstrated in the inset in the right frame of 
Fig. 8. Notice that these oscillating potentials will also have pronounced effects 
on the interaction between such rigid polyclcctrolytes. HNC/MSA theory can 
qualitatively and quantitatively describe such effects to a very good degree, as 
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has been demonstrated theoretically a long time ago [32]. A recent comparison 
with simulations can be found in Ref . [36] . 



5 Pressure 

5.1 Pressure measurements within the generic model 

Using the results from the Appendix A we now turn towards pressure measure- 
ments for the generic systems described above. Because of the periodic bound- 
ary conditions, we are actually looking at the pressure of a nematic, hexagonally 
packed solution of charged rods. We specifically compute the dimensionless os- 
motic coefficient p, which is simply the pressure normalized by its ideal gas 
contribution 

^ P (16) 
p's nk^T 

For an isolated cell the pressure is given by the particle density at the outer 
cell boundary. This is a rigorous statement, true for the spherical, cylindrical 
and planar cell model [9]. It is a merit of the PB equation that it retains the 
validity of this exact relation. For the extended density functional theories this 
no longer holds and additional terms appear [29] . 

Within Poisson-Boltzmann theory the osmotic coefficient as determined from 
the boundary density is 

where 7 is the density dependent integration constant from Eqn. (11), 

Fig. 9 illustrates these measurements for monovalent counterions and three 
values of ^b/c = 1,2,3. Several things may be noted: The osmotic coefficient 
from the simulations is always smaller than the PB prediction, but for low den- 
sity both values converge. This also illustrates that the Manning limiting law 
from the r.h.s of Eqn. (17) becomes asymptotically correct for dilute systems. 
Upon increasing the density, the osmotic coefRcient rises weaker than the PB 
prediction. This is more pronounced for systems with higher Bjerrum length, 
and consequently, higher Manning parameter, and is due to enhanced coun- 
terion condensation which has been observed in Sec 4. Notice that this has a 
very remarkable side-effect. Over a considerable range of densities the measured 
osmotic coefficient is much closer to the limiting law than to the actual PB pre- 
diction. This makes the Manning limit look much more accurate than it really 
ought to be. However, the surprising effect should not be over-interpreted, since 
the underlying reason is nothing but a fortunate cancellation of two contribu- 
tions of approximately the same size which are not contained in the limiting 
laws. 

Also for the pressure it is interesting to investigate complementary systems in 
which the values of Bjerrum length £b/(J and valence v have been interchanged 
to keep the product fixed. However, at given density the cell radius depends 
on the valence, so p{n) does no longer universally depend on this product. Since 
in the limit of high density p is easily seen to increase with density, while at low 
density it decreases, the functions p{n) for different valences intersect at some 
point. 
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Fig. 10 summarizes the results of measurements on the multivalent systems 
with V = 1,2, 3, which yield the same values of as the monovalent ones inves- 
tigated before. The most striking feature is the appearance of a negative osmotic 
coefficient in a certain density region of the trivalent case. If the constraint of 
a fixed rod-separation were to be replaced, the system would phase separate, 
hence, attractive interactions must be present between the rods. Similar ob- 
servations have been reported in Rcfs. [24, 25, 26]. In order to trace back the 
origin of those attractions, Fig. 11 displays the osmotic coefficient in the diva- 
lent and trivalent cases, split into two contributions: (i) the non-electrostatic 
part comprising ideal gas contribution and the short-range virial and (ii) the 
electrostatic part, which for visual convenience is plotted with a reversed sign. 
Hence, the difference between those two curves gives the curves in Fig. 10. We 
observe that the electrostatic part leads to a monotonically increasing attraction 
with increasing density. It is almost twice as strong in the trivalent case. The 
very strong increase of the osmotic coefficient at large densities is due to the 
virial, i.e., due to repulsive ion-rod and ion- ion interactions. The measured neg- 
ative pressure in the trivalent case is the result of a "sudden" drop in the virial 
contribution. Nothing particular can be observed in the electrostatic part. We 
suggest the explanation that in an intermediate range of densities the presence 
of surrounding rods lifts condensed ions from the surface they are pushing on, 
thereby reducing the short range excluded volume contribution. In too dilute 
systems this effect is not operative, while in too dense systems the effect of the 
other rods is actually to push ions more closely to the surface. Hence, there is a 
density range in which the virial is reduced by the presence of other rods. 

Contrary to the simulations, the osmotic coefficient from the PB theory is 
always positive. This is the consequence of the rigorous proof that such attractive 
interactions are absent on the PB level [37]. An extension of this statement 
to ions of finite size and to a wider class of boundary conditions and density 
functional can be found in Ref. [38]. 

Finally it should be noted that the above measurements can not be used to 
infer that attractive forces between charged rods require the counterions to be 
at least trivalent. The reason is twofold: First, at given valence one can vary 
Bjerrum length and line charge density. Increasing the Manning parameter will 
lead to negative pressure in the divalent system. Second, keeping all interaction 
potentials fixed, the radius tq of the charged rod is a relevant observable, as has 
been demonstrated in Ref. [25]. Hence, a general statement about presence or 
absence of attractive interactions is difficult, since a five-dimensional parameter 
space is involved: {X,iB,v,ro,n}. 

6 Examples: poly(p-phenylene) 

This section discusses simulations in which the parameters are explicitly mapped 
to an experimental system. This affects the values of ion diameter a, rod radius 
rp, Bjerrum length £b and line charge density A. In order to have rg different from 
(T, we use as the ion-rod-potential a LJ like expression, in which r is replaced 
by r — Tg and the shift rg is related to the desired rod radius by ro = rg 4- cr/2. 

Since these systems have been investigated either in the presence of salt or 
at rather high Manning parameter, the phenomena happening in the vicinity of 
the rod are always well decoupled from the cell boundary. Therefore, an even 



11 



simpler geometry has been chosen: One rod placed parallel to one of the edges 
of a cubic simulation cell. 

6.1 Distribution functions and osmotic coefficient for salt 
free poly(p-phenylene) solutions 

The Manning limiting laws are claimed to apply in the limit of low density. 
However, an experimental verification of this effect using DNA as rod-like poly- 
electrolytes is difficult, since the double helix starts to unwinds at low ionic 
strength. It would therefore be desirable to have a stiff model polyelectrolyte 
which docs not suffer from this problem. One such system belonging to the class 
of poly(para-phcnylcnes) (ppp) has recently been investigated in Refs. [16, 17]. 
The constitution formula is given in Fig. 1. The fully aromatic backbone ex- 
hibits an excellent chemical stability and has a persistence length of approxi- 
mately 20 nm. The degree of polymerization used in the abovementioned studies 
was between 20 and 40, so that the contour length equals approximately one 
persistence length at most. The dominant contribution to the signal in small 
angle X-ray scattering experiments stems from the iodine ions, since the excess 
electron density of the backbone is very low. Tabic 1 lists four systems which 
have been simulated based on a mapping to ppp. 

In a first step the simulated ion distribution functions shall be compared 
with the PB prediction as well as with the Debye-Hiickel-hole-cavity theory from 
Ref. [15], which tries to incorporate correlations missing in the PB approach. 
Figure 12 shows the corresponding distribution functions for the systems from 
Table 1. Qualitatively the PB prediction is already a good description of those 
systems, if one graciously disregards the inevitable problems at the cell bound- 
ary. However, the presence of correlations shifts the distribution functions up 
in all four cases. Since this shift is rather small, i.e., the correlations are weak, 
the Debye-Hiickel-hole-cavity theory is in fact an excellent description of those 
systems. Observe, for instance, that in the weakly charged dilute case its pre- 
diction can no longer be distinguished from the simulated curve on the chosen 
scale. 

With the help of the condensation criterion from Sec. 3.2 the extent of the 
correlation-enhanced condensation can be further quantified. The Manning frac- 
tion from the molecular dynamics simulation increases only by a fairly small 
amount. It is at most 4% larger than the PB value. This translates to an ef- 
fective Manning parameter ^eff = 1/(1 ^ /m) being 5-10% larger than the bare 
one. This increase is very accurately captured by the Debye-Hiickel-hole-cavity 
theory. Its prediction for /m is at most 1% smaller than the value obtained in the 
simulation. Interestingly, it is even independent of density. This, however, is not 
a feature to be generally expected and should therefore not be over-interpreted. 

In a second step the osmotic coefficients of the four model systems have been 
computed. Two distinct approaches have been used for analyzing the simula- 
tion results: The first computes the pressure from the stress tensor, as described 
in detail in Appendix A. The second approach exploits the connection from 
Eqn. (17) between osmotic coefficient. Manning parameter and the integration 
constant 7 of the PB equation. The idea is to use the fitting procedure described 
in Sec. 4.2 and from that obtain values and 7*, leading to the osmotic co- 
efficient (1 + 7*^)/2^*. The most straightforward approach of measuring the 
boundary density is less recommcndablc. since in the simulation the boundary 
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has a quadratic and not a circular cross section, i.e., the merely approximate 
representation of the cell model becomes most visible and questionable here. In 
contrast to that, the two other approaches "look" at regions of the cell, which 
are away from the boundary. For comparison, also the predictions from the 
Manning limiting law and the full PB expression are calculated. For the two 
dilute systems there also exist measurements in Ref. [17] of p via osmometry, 
while for the dense systems this was unfortunately impossible due to counter 
diffusion problems. Tabic 2 collects the coefficients for comparison. 

The first thing which should be noted is that PB theory gives again a surpris- 
ingly good description of the measured coefficients, including the experimentally 
determined one. A closer look reveals that - as expected - it overestimates p, 
and the deviations are larger for the dense systems. Still, PB theory is off by 
only 7% or 18% for the dilute or dense systems, respectively, thereby confirm- 
ing the observation already made when looking at the distribution functions. A 
further indication of this point is that the Manning limit 1/2^ appears to be a 
lower boundary, which is also in accord with the pressure measurements on the 
generic systems with monovalent ions from Sec. 5.1. Since there it has also been 
found that the Manning limit need no longer act as a boundary in the multi- 
valent case, see Fig. 10, it would be very interesting to perform experiments on 
those systems with, e.g., divalent ions. Observe that the experimental value of 
p for the highly charged dilute system 2 is already at the Manning limit. 

Since the experimentally determined osmotic coefficient appears to be smaller 
even than the molecular dynamics results, this indicates effects to be relevant 
which go beyond the model used for simulation. Most obvious candidates for 
this are the neglect of additional chemical interactions between the ions and the 
polyelectrolyte as well as solvation effects, i.e., interactions between the ions or 
the polyelectrolyte with the water molecules from the solution. It is for instance 
demonstrated in Ref. [17] that the osmotic coefficient also depends on whether 
on uses chlorine or iodine counterions. While one could certainly account for the 
different radii of these ions when computing the distance of closest approach en- 
tering the PB equation, the implications of the different hydration energies is 
much less obvious to incorporate and in principle requires very expensive all- 
atom simulations. See also the discussion in Ref. [18], which includes the effects 
of small excess salt to the osmotic coefficient. 

Finally it should be noted that the two different methods of analyzing the 
molecular dynamics data lead to very similar results. While this fact is not par- 
ticularly useful in a simulation, its consequences for the analysis of experimental 
data seem promising. In Ref. [16] poly(p-phenylene) solutions are investigated 
by means of small angle X-ray scattering, which due to the rod-like geometry 
turns out to be sensitive to the radial ion distribution function. Therefore, the 
measured structure factors can be related to PB distribution functions with ef- 
fective values of ^, 7 etc. Fitting the scattering intensity and thereby determining 
those values permits in principle a measurement of the osmotic coefficient along 
the lines already described above. What makes this approach so attractive is 
that it could ideally complement osmometry. If the density becomes too large, 
the latter suffers from severe problems with counter diffusion, but it is exactly 
this high density which meets the requirements of good contrast necessary for 
scattering. 

A similar fitting procedure has been applied to the measured data from small 
angle X-ray scattering experiments on the systems 3 and 4; see also Ref. [16]. 
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From a fit to the structure factor the radius Tq of the rod has been determined. 
The integration constant 7* and the Manning parameter i?^ the follow from 
the PB theory. The obtained values are listed in Table 3. The resulting Manning 
radius is surprisingly close to the one determined by simulation, particularly for 
the system with ^ = 3.32. The osmotic coefficients constructed from the PB 
formula are less accurate, 10-20% too large. 

7 The importance of correlations 

We have seen that the nonlinear PB equation suffers from systematic deviations 
in strongly coupled or dense systems. It underestimates the extent of counterion 
condensation and at the same time overestimates the osmotic coefficient. The 
common reason for both problems is the neglect of correlations: More elaborate 
theories [32, 33, 15] which try to incorporate correlations reproduce the trends of 
the simulations, because they favor an increase in density close to the macroion 
surface, thereby leading to a stronger condensation and a concomitant drop in 
the boundary density and thus pressure. 

On the one hand there is general agreement about the presence of correla- 
tions being liable for the failure of PB theory. On the other hand, this insight 
does not shed any light onto the question, which kind of correlations are impor- 
tant. In the following we investigate a simple pair-correlation function based on 
a one-rod property. 

If the surface charge density of the rod is high, a fairly large number of 
countcrions will stay within a condensed layer of small radial extent. Addition 
of salt will further increase the ionic density in this layer. It is clear that beyond 
a certain point (high A, high £b or enough salt) the ions will no longer distribute 
independently of each other but get locally correlated. This effect will now be 
measured for a DNA-like system, because for those systems correlation effects 
are assumed to play a prominent role. The rod radius was chosen to be 7.86 
A(1.85(t) and the ion diameter is correspondingly 4.25 A(1ct). In addition 0.5 
mol/1 of a 2:2 salt has been added in excess to the divalent countcrions. More 
information about theses systems can be found in Table 4. 

The first issue is to define the innermost layer. For this, a distance from the 
rod is chosen which contains many countcrions but virtually no co-ions. A dis- 
tance of roughly 11.5 A from the rod axis turned out to be suitable. This is about 
a third ion diameter farther out than the distance of closest approach. To avoid 
difficulties with remaining co-ions, only the countcrions within this distance are 
taken into account in all what follows. In a second step the coordinates of those 
ions are radially projected onto the surface of the cylinder of closest approach, 
and this surface is then rolled out to a flat plane, see Fig. 13 for an illustration 
of this procedure. Finally, the two-dimensional pair correlation function g{r) of 
these projected points is computed. 

One might object that the unrolled flat plane leads one to believe that the 
pair correlation function does only depend on distance, while in reality such a 
rotational symmetry cannot be expected, due to the rod curvature. However, 
an investigation of g{r) revealed that the pair correlation function is in fact 
rotationally symmetric up to essentially a distance which corresponds to half 
the circumference of the cylinder of closest approach, which is roug hly 30 A. 
Larger distances can of course only be realized along the rod instead of around 
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it. A possible reason for this surprising symmetry is that at short distances the 
curvature is not yet perceptible while at large distances the particles are already 
uncorrelated. 

Figure 14 shows the measured g{r) for four differently charged systems with 
Manning parameter 10.5, 8.4, 6.3, and 4.2. The last value corresponds to DNA 
in aqueous solution. The most important thing to observe is that for the strongly 
charged systems g{r) shows definite signs of correlations. Apart from the trivial 
correlation hole at small r there is a distance rmax at which ions are more likely 
to be found than under the assumption of independent distribution. Notice that 
the maximum in g{r) is not located at r « cr but much farther out. Hence, its 
existence is not merely an artifact of close packing of repulsive Lcnnard- Jones 
spheres. However, if the condensed ions are assumed to form a triangular lattice 
on the surface in order to maximize their mutual repulsion, the resulting distance 

is 25 - 35 % larger than rmax- Together with the only weakly pronounced 
oscillations in g{r) this proves the correlation induced interactions to be rather 
short-ranged, yet less local than a pure hard core. In any case, the range is 
several times larger than the average salt Debye length £d w 0.5 ct. 

Such two-dimensional pair correlation functions for mobile ions adsorbed on 
charged planes have recently been investigated theoretically in Ref. [39]. The 
mean separation = {veo/<;y^'^ of w-valent ions on the completely neutralized 
plane of surface charge density is identified as the important scaling length. 
The main predictions for the strongly coupled regime r2 = t^v^ /dc^ > 1 are: 
(i), g{r) should have a single first peak at a distance about the size d^ of the 
correlation hole, (ii), the breadth of this peak should decrease with r2 while its 
maximum should increase. This trend is indeed seen to be true for the functions 
in Fig. 14, only d^ is roughly 10-20% larger than the actual peak position. One 
might want to argue that the presence of much salt entails a further increase 
of the layer density on top of the usual Gouy-Chapman prediction, but the 
results obtained in Ref. [39] are claimed to be independent of the bulk ionic 
strength, if the latter is smaller than the layer density. This is here the case, even 
though not always by an order of magnitude. An alternative explanation may be 
based on the presence of co-ions: Although there reside very few of them within 
the innermost condensed layer, there will be an appreciable amount beyond 
and possibly very close to n- Those ions can act as "bridges", they attract 
counterions and thereby reduce the average closest distance between them. 

Although Tmax < ^ A , this does not exclude a local hexatic ordering of the 
counterions. But since g{r) is on average rotationally symmetric, establishing 
signs for this requires the investigation of a suitable 3-point correlation func- 
tion. The trick is to break the rotational symmetry, which suggests the following 
procedure: g(r) is proportional to the probability of finding a particle at a dis- 
tance r, given that there is also a particle at the origin. Now define g^{r) to 
be proportional to the probability of finding a particle at position r, given that 
there is also a particle at the origin and given that the particle which is closest 
to the origin is situated to the right side. This definition is further restricted to 
be sensitive only for next nearest neighbors and not arbitrary other particles. 
Thence, it answers the question: "If the nearest neighbor is on the right side, 
where is the next nearest neighbor?" A scatter plot of this observable is shown 
in the left part of Fig. 15. For this, the most highly charged DNA-sized system 
with Manning parameter ^ = 10.5 has been used. Clearly visible are the two cor- 
relation holes around the origin and the nearest neighbor as well as a tendency 
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of the next nearest neighbor to be on the side opposing the nearest neighbor. 
All those effects are trivially explained. However, on top of that no preferential 
ion accumulations in the "hexatic directions" indicated as dotted lines is per- 
ceptible. The right part of Fig. 15 shows the probability distribution of finding 
the next nearest neighbor at an angle cj) with respect to the line joining nearest 
neighbor and origin. The correlation hole generated by the nearest neighbor is 
visible, but no increase in p{(j)) at 60°, 120° or 180°. Rather, the measured data 
are compatible with a fairly structureless distribution, as indicated by the solid 
line. 

Since the most strongly charged system does not show hexatic order, it will 
certainly be absent in the other ones. However, the coupling constant T2 intro- 
duced above has a value only about 3 in the system with Manning parameter 
^ = 10.5. For larger values of T2 one would not only expect much stronger corre- 
lations and hexatic ordering but even crystallization of the adsorbed counterions 
into a two-dimensional Wigner crystal [39, 40]. More extensive investigations to 
clarify the role of correlations are presently under way [41] 

8 Summary 

Theoretical and numerical studies of stiff linear polyelectrolytes within the 
framework of a cell model have been presented. We outlined methods to quan- 
tify condensation for general measurements of ion distribution functions via an 
inflection point criterion and by using a generalized PB fit function. We com- 
pared our simulational results to predictions from mean-field PB theory, and 
found excellent agreement for weakly charged systems. The enhanced counter- 
ion condensation caused by ionic correlations and its dependence on parameters 
such as density, Bjerrum length, valence and ionic strength was measured, and 
effects which qualitatively go beyond mean- field theory, e.g., charge reversal and 
attractive interactions between like-charged macroions have been found. For sys- 
tems with a high salt content, integral equation theories have proved to give a 
good qualitative and quantitative description of such systems [32, 36] 

For the osmotic pressure measurements we displayed for a large density range 
that the measured pressure stays close to the Manning limiting law prediction 
due to a fortunate error cancellation of neglected ionic correlations and density 
effects which point in different directions. We showed for the specific example 
of a ppp-polyelectrolyte what differences to a PB description can be expected 
for both, the ionic distribution and the osmotic pressure. We finally analyzed 
some correlations for a DNA-like system, and showed that they show signs of a 
strongly correlated hquid, as has been advocated recently [39, 40]. 

To assess these theoretical ideas, computer simulations will certainly be in- 
dispensable, since they are presently the only practicable way of obtaining suf- 
ficiently detailed information on ionic distributions and correlation functions. 
The numerical investigations presented in this work and the correlation analysis 
of the obtained data are a further step in this direction, however a more detailed 
analysis is presently under way [41] . 
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Appendix: 

A Defining and computing the pressure 

The pressure for the primitive cell model is a nontrivial variable to compute for 
two reasons: First, the long-range electrostatics has to be properly taken into 
account, and secondly, the system is inherently anisotropic. Hence, the relevant 
observable is the stress tensor. 

For isotropic systems, the thermodynamic definition of the pressure as the 
derivative of the free energy with respect to volume leads to the general equation 

pV = NkBT-v(^^y (18) 

where the angular brackets (• • • ) denote a canonical average. For the case of 
short-range interactions the contribution from U can be further simplified by 
using 

dU v-^ dU dri v-^, „ , , , 

ay =Ea;r-^-E(-^0-^. (19) 

Substituting this into the pressure equation (18) and using Fij ~ ^^ji gives 

- iVfeeT-i ^(ry -F,,) (20) 

i<j 

The electrostatic contribution to the pressure can be derived [42] by realizing 
that the electrostatic energy is a homogeneous function of volume with 
degree —1/3, and hence by Euler's theorem dU'^/dV = —U'^/3V. Adding this 
to Eqn. (20) finally gives the desired pressure equation 

pV = 7VfcBT-i^(r,,.F-) + i (C/C), (21) 

i<j 

where F^j are the short-range pair forces. Note that for a repulsive hard core the 
virial contribution is positive, while the electrostatic contribution is one third 
of the energy density, which for neutral systems is usually negative. 

For anisotropic systems the stress tensor is the relevant observable to com- 
pute. Whereas the ideal gas contribution to the pressure still remains isotropic, 
the scalar product in the virial is replaced by the tensor product "(g)" . For the 
case of electrostatic interactions the derivation is more complicated and will not 
be presented here. The reader is referred to Ref. [43]. The result is that within 
the framework of Ewald techniques the electrostatic contribution to the stress 
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tensor, p*-', can be decomposed additively into a real space contribution p*^*"^ and 
a Fourier space contribution p*-*^^. They are given by 




(22) 



(23) 



Here, rijm ~ fi ~ rj + mL\^, and the canonical average has not been denoted 
explicitly. For details of the notation and the Ewald techniques in general, the 
reader is referred to the literature [21, 22]. Note that the connection to the 
isotropic case requires the trace of the pressure tensor to be equal to the mean 
electrostatic energy density, i.e., Trp^ = U'^/V. 

After computing the electrostatic stress tensor, it can be rotated such that 
the rod points along the z-axis. The xx- and yy-components of the new tensor 
then give the pressure perpendicular to the rod, while the zz-component is the 
contribution parallel to the rod. 
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List of Tables: 



system 


c [g/1] 


T [°C] 




a [A] 


rfca [A] 


R[k] 


1 


1.5 


40 


3.4 


7.9 


4.4 


239.3 


2 


1.5 


40 


6.8 


8.1 


4.4 


292.0 


3 


19.95 


25 


3.32 


7.9 


4.4 


65.6 


4 


17.99 


25 


6.64 


8.1 


4.4 


84.3 



Tabic 1: Mapping for four poly(p-phcnylene) systems. Tabulated are polyelec- 
trolyte concentration c, temperature T, Manning parameter ^, ion diameter cr, 
distance of closest approach c?ca between ions and the rod and cell radius R cor- 
responding to the given concentration. The experiments have been performed 
under salt-free aqueous conditions with monovalent countcrions [16]. 



system 




Ppb 


Pmd,i 


PMD,2 


P cxp 


1 


0.1471 


0.2128 


0.2005(57) 


0.201 


0.185(15) 


2 


0.0735 


0.1073 


0.1003(60) 


0.102 


0.073(15) 


3 


0.1506 


0.2848 


0.2424(57) 


0.260 




4 


0.0753 


0.1428 


0.1215(56) 


0.129 





Table 2: Osmotic coefficient of the four poly(p-phenylene) systems from Table 1. 
Shown are the Manning limiting law (M). the PB prediction, the simulation 
result using the stress tensor (MD,1), the simulation result using a PB fit as 
described in Sec. 4.2 (MD,2) and results from measurements on this system 
using osmometry [17] 
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system 




R[k] 


^0* [A] 


R*M [A] 


7* 


P 


3 


3.32 


65.6 


8.2 


28.19 


0.957 


0.289 


4 


6.64 


84.3 


9.9 


39.14 


1.014 


0.153 



Table 3: Fit evaluation of small angle X-ray scattering experiments on the sys- 
tems 3 and 4 [16]. ^ is the Manning parameter, R the cell radius implied by the 
polyelcctrolyte concentration, and Tq the rod radius, which had been the only 
fitting parameter. From this, the integration constant 7* and the Manning ra- 
dius is determined from the solution of the PB equation, p* = (1 -I- 7*^)/2^ 
is the attempt to predict the osmotic coefficient of the solution with the help of 
the PB formula. 



parameter 


symbol 


value 


value in LJ units 


ion diameter 


a 


4.25 A 




ion valence 


V 


2 


2 


rod radius 


ro 


7.86 A 


1.85 a 


line charge density (DNA) 


X 


eo/l.7A 


2.5 eo/cr 


Bjerrum length (water) 


(b 


7.14 A 


1.68 a 


Manning parameter 


c 


4.2 


4.2 


box size 


Lb 


122.4 A 


28.8(7 


corresponding cell radius 


R 


69.1 A 


16.2 a 


temperature 


T 


298 K 


e/k-B 



Table 4: Simulation parameters for a DNA-likc system. 
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List of Figures 




Figure 1: Example of a stiff polyelectrolyte. Constitution formula for poly(para- 
phenylene) with iodine counterions (left) and a physicist's picture (right). 
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Figure 2: Realization of the cell model. A rod placed along the main diagonal 
of a cube yields an infinite triangular array of infinitely long rods upon peri- 
odic replication of the original cube. The Wigncr-Scitz cell of this lattice is a 
regular hexagon enclosing the rod. This therefore provides a way of obtaining a 
hexagonal cell without abandoning a cubic geometry. 




Figure 3: Countcrion distribution functions -P(r), as defined in Eqn. (13), for 
various simulated densities. Note that an increasing cell radius corresponds to 
functions extending towards larger values of r. The heavy dots mark the points 
of inflection in P as a function of ln(r), while the crosses mark the positions 
at which those points would be located on the corresponding PB distribution 
functions. For the sake of clarity such a PB distribution is only plotted for the 
system with lowest density, i.e., R = 123.85(7 (dotted line). 
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Figure 4: Counterion distribution functions P{r) (solid lines) for seven systems 
with the same dimensions as the ones in Fig. 3, but with a Bjcrrum length 
^b/c = 1- Since the resulting Manning parameter ^ = 0.959 < 1, counterion 
condensation is not expected to occur. This is borne out by the observation 
that the functions are convex up already at r = rg. In these weakly charged 
systems the predictions of PB theory (dotted lines) are excellent and can hardly 
be distinguished from the simulation results. 
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Figure 5: The functional form in Eqn. (14) has been fitted to the measured 
distribution hmction (solid line) of the system with R = 123.85 ct and £b = 
3(7 within the range [4ro;40ro] between the two vertical bars. The |-arrow 
indicates the inflection point of the PB distribution (dotted line), while the [- 
arrow indicates the corresponding point (-Rmj/m) ~ (11.0 cr; 0.664) of the fit 
(gray stripe). Note that i?^ < i?M although > /m- The result of the fit is 
C = 2.97, 7* = 0.457 and = 0.579 ct. 




Figure 6: Counterion distribution functions P{r) versus r/ro. The high (low) 
density situation is shown in the left (right) frame. The lower three curves are 
for vi^ja = 2, while the upper three correspond to vI-q/g = 3. The systems 
with multivalent counterions (solid lines) always show a stronger condensation 
than the complementary systems with monovalent ions (dashed lines), which 
themselves show a stronger condensation than PB theory (dotted lines) . 
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Figure 7: Distribution functions P{r) for a system which has R/<t = 61.923, 
^b/c = 2.189, ^ = 2.1 and monovalent ions. From bottom to top the number 
of sah molecules added to the simulation box of length Lb = 144.446 a is 0, 
104 and 3070, which corresponds to a salt Debye length of oo, 22.9(7 and 4.2 cr, 
respectively. The solid lines are the result of a simulation while the dotted lines 
are the predictions of PB theory. 
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Figure 8: Charge distribution for a system characterized by Lh/a = 36.112, i.e., 
60 monovalent counterions and a cell size of R/a = 15.481, l^ja = 4.1698, i.e., 
^ = 4, and 1000 molecules of a 2:2 salt within the box volume. The simulation 
shows a pronounced overcharging-effect in P(r) (solid curve, left frame), in 
contrast to PB-theory (dotted curve). The charge oscillations can be described 
quite accurately by an exponentially damped sine function with period 1.89 cr 
and decay length 0.85 cr (gray stripe). The densities n-2{r) (solid line, right 
frame) and n_|_2(r) (dotted line) of negative and positive salt ions, respectively, 
demonstrate the effect of charge layering and local charge reversal, and the 
inlay shows the dimensionless electrostatic potential y(r) — (3eQip{r), which is 
also oscillating. 
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Figure 9: Osmotic coefficient p versus na^ for monovalent countcrions. Heavy 
dots mark the measurements, while the solid lines are fits which merely serve 
to guide the eye. The dotted lines are the prediction of PB theory. From top to 
bottom the Bjerrum length Ib/ct varies as 1,2,3. The errors in the measurement 
are roughly as big as the dot size. 
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Figure 10: Osmotic coefficient p as a function of density n for different valences. 
Heavy dots mark tlic measurements, while the solid lines arc fits which merely 
serve to guide the eye. The dotted lines arc the prediction of PB theory. From 
top to bottom the counter ion valence v varies like 1,2,3, which gives the same 
value of as the curves in Fig. 9. The errors in the measurement are roughly 
as big as the dot size. 
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Figure 11: Osmotic coefficient p for the divalent (left) and trivalent (right) sys- 
tems from Fig. 10, separated into the non-electrostatic contribution coming from 
virial and ideal gas (heavy dots on solid lines) and negative electrostatic contri- 
bution (crosses on dotted lines). Again, the lines are fits which merely serve to 
guide the eye. 
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Figure 12: Distribution functions for the 4 poly(p-phenylene) systems from Ta- 
ble 1. Left/right frame eorrespond to the high/low density systems (3,4)/(l,2), 
while the upper/lower set of functions correspond to the strongly/weakly 
charged systems (2,4)/(l,3) respectively. Solid lines are the results of simula- 
tion, dotted lines are the PB prediction and dashed lines are from an extended 
PB theory using the Debye-Hiickel-hole-cavity correction [15]. The "f-arrows in- 
dicate the inflection points in the PB distributions, while the j-arrows mark 
those points in the simulated distributions. The deviations of the latter from 
the PB curves at large r originate from the simulation cell having a quadratic 
instead of a circular cross section. 




Figure 13: Illustration of the computation of surface correlations. Counterions 
within a certain small distance from the rod constitute the innermost condensed 
layer. Their coordinates are radially projected onto the surface, which after that 
is unrolled to a flat plane. The two-dimensional pair correlation function g{r) is 
then determined from the projected points. 
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Figure 14: Two-dimensional pair correlation function g(r) for the projection of 
counterions within a close condensed layer onto the cylinder of closest approach, 
see text and Fig. 13. For the four functions the Manning parameter decreases 
as 10.5, 8.4. 6.3 and 4.2. The last value corresponds to DNA. 
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Figure 15: Left part: Density plot of the 3-point correlation function (^(r) for 
the highest charged system from Fig. 14 with Manning parameter ^ = 10.5. If 
an ion is located at the origin and its nearest neighbor is on the right side and 
not further away than 1.3 cr, its next nearest neighbor is shown as a dot. The 
circle around the origin has radius a and corresponds to the distance of closest 
approach imposed by the Lennard- Jones potential. The right frame plots the 
probability distribution of finding the next nearest neighbor at an angle </> with 
respect to a line joining nearest neighbor and origin. The solid line is a guide to 
the eye. 
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